#include "mbpt.h"

vector<int> operator+(vector<int> &v1, vector<int> &v2) {
  // 假设 v1.size() == v2.size()
  if (v1.size() != v2.size()) {
    cout << "Wrong @ vector<int> operator+ " << endl;
    exit(0);
  }
  vector<int> r;

  r.reserve(v1.size());
  for (auto i = 0; i < v1.size(); ++i) {
    r.push_back(v1[i] + v2[i]);
  }
  return r;
}
vector<int> operator+(const vector<int> &v1, const vector<int> &v2) {
  // 假设 v1.size() == v2.size()
  if (v1.size() != v2.size()) {
    cout << "Wrong @ vector<int> operator+ " << endl;
    exit(0);
  }
  vector<int> r;

  r.reserve(v1.size());
  for (auto i = 0; i < v1.size(); ++i) {
    r.push_back(v1[i] + v2[i]);
  }
  return r;
}

vector<double> operator+(vector<double> &v1, vector<double> &v2) {
  // 假设 v1.size() == v2.size()
  if (v1.size() != v2.size()) {
    cout << "Wrong @ vector<double> operator+ " << endl;
    exit(0);
  }
  vector<double> r;

  r.reserve(v1.size());
  for (auto i = 0; i < v1.size(); ++i) {
    r.push_back(v1[i] + v2[i]);
  }
  return r;
}

vector<double> operator+(const vector<double> &v1, const vector<double> &v2) {
  // 假设 v1.size() == v2.size()
  if (v1.size() != v2.size()) {
    cout << "Wrong @ vector<double> operator+ " << endl;
    exit(0);
  }
  vector<double> r;

  r.reserve(v1.size());
  for (auto i = 0; i < v1.size(); ++i) {
    r.push_back(v1[i] + v2[i]);
  }
  return r;
}

vector<double> operator-(vector<double> &v1, vector<double> &v2) {
  // 假设 v1.size() == v2.size()
  if (v1.size() != v2.size()) {
    cout << "Wrong @ vector<double> operator+ " << endl;
    exit(0);
  }
  vector<double> r;

  r.reserve(v1.size());
  for (auto i = 0; i < v1.size(); ++i) {
    r.push_back(v1[i] - v2[i]);
  }
  return r;
}

void MBPT_SEC::init() {
  this->sp_core = norm_ct.sp_core;
  this->sp_core_Tz = norm_ct.sp_core_Tz;
  this->sp_P_beg = norm_ct.sp_P_beg;
  this->sp_P_beg_Tz = norm_ct.sp_P_beg_Tz;
  this->sp_P_end = norm_ct.sp_P_end;
  this->sp_P_end_Tz = norm_ct.sp_P_end_Tz;
  this->sp_Q_beg = norm_ct.sp_Q_beg;
  this->sp_Q_beg_Tz = norm_ct.sp_Q_beg_Tz;
  this->sp_Q_end = norm_ct.sp_Q_end;
  this->sp_Q_end_Tz = norm_ct.sp_Q_end_Tz;
  this->hbar_omega = norm_ct.hbar_omega;
  //this->E_cut = norm_ct.E_cut;

  Frac_vec.resize(Derivative_Max + 1);
  Frac_vec[0] = 1;
  for (int i = 1; i <= Derivative_Max; i++) {
    Frac_vec[i] = Frac_vec[i - 1] * i;
  }

  cout << " ---- MBPT_SEC Parameter ---- " << endl;
  cout << "\t hbar_omega = " << hbar_omega << endl;
  cout << "\t E_cut = " << E_cut << endl;

  cout << "\t sp_core_Tz = " << sp_core_Tz << endl;
  cout << "\t sp_P_beg_Tz = " << sp_P_beg_Tz << endl;
  cout << "\t sp_P_end_Tz = " << sp_P_end_Tz << endl;
  cout << "\t sp_Q_beg_Tz = " << sp_Q_beg_Tz << endl;
  cout << "\t sp_Q_end_Tz = " << sp_Q_end_Tz << endl;

  this->Energy_shift = norm_ct.Energy_shift;

  this->sixJ_x.init(12);
  this->Sqrt_vec.resize(50);
  for (int i = 1; i < 50; i++) {
    double t = i;
    Sqrt_vec[i] = std::sqrt(t);
  }
  cout << "Start to build_Tz_Half" << endl;
  DP_Model.build_Tz_Half(SP, sp_P_beg_Tz, sp_P_end_Tz, sp_P_beg_Tz,
                         sp_P_end_Tz);
  cout << "Start to print_Tz" << endl;
  DP_Model.print_Tz("data/MBPT_DP.dat");
}

void MBPT_SEC::build_sec_Sm() {
  cout << " this->hbar_omega = " << this->hbar_omega << endl;
  int sp_size = SP.state_Tz.size();
  Sm_vec.resize(sp_size);

  for (int i = 0; i < sp_size; i++) {
    Sm_vec[i].resize(Derivative_Max + 1);
  }
  for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {

    vec_1B[a] = cal_OneBody(a, Sm_vec[a]);
    double val_fist = get1B(a, a);
    Sm_vec[a][0] += val_fist;
    for (int i = 0; i <= Derivative_Max; i++) {
      // Times m!
      Sm_vec[a][i] = Sm_vec[a][i] * Frac_vec[i];
    }
  }
  this->Sm_build_flag = 1;
}

void MBPT_SEC::print_sec_OB_sep() {
  cout << " this->hbar_omega = " << this->hbar_omega << endl;
  int sp_size = SP.state_Tz.size();
  vec_OA.resize(sp_size);
  vec_OB.resize(sp_size);
  vec_OC.resize(sp_size);
  vec_OD.resize(sp_size);
  vec_OE.resize(sp_size);
  vec_1B.resize(sp_size);
  for (int i = 0; i < sp_size; i++) {
    vec_OA[i] = 0.0;
    vec_OB[i] = 0.0;
    vec_OC[i] = 0.0;
    vec_OD[i] = 0.0;
    vec_OE[i] = 0.0;
    vec_1B[i] = 0.0;
  }
  vector<vector<vector<double>>> vec_d;
  vec_d.resize(sp_size);
  for (int i = 0; i < sp_size; i++) {
    vec_d[i].resize(5);
    for (int k = 0; k < 5; k++) {
      vec_d[i][k].resize(Derivative_Max + 1);
      for (int j = 0; j <= Derivative_Max; j++)
        vec_d[i][k][j] = 0.0;
    }
  }
  for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
    // cout << "\t a = " << a << endl;
    double E_start = SP.state_Tz[a].E * hbar_omega;
    vec_OA[a] = cal_OA(a, E_start, vec_d[a][0]);
    vec_OB[a] = cal_OB(a, E_start, vec_d[a][1]);
    vec_OC[a] = cal_OC(a, E_start, vec_d[a][2]);
    vec_OD[a] = cal_OD(a, E_start, vec_d[a][3]);
    vec_OE[a] = cal_OE(a, E_start, vec_d[a][4]);
    vec_1B[a] = vec_OA[a] + vec_OB[a] + vec_OC[a] + vec_OD[a] + vec_OE[a];
  }
  // cout << "a\t A\t\t B\t\t C\t\t D\t\t  E*2\t " << endl;
  for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
    cout << " --- SP Tz : " << a << "\t ---" << endl;
    cout << "\t A : ";
    for (int i = 0; i <= Derivative_Max; i++)
      cout << std::setw(8) << std::setprecision(5) << vec_d[a][0][i] << "\t ";
    cout << endl;
    cout << "\t B : ";
    for (int i = 0; i <= Derivative_Max; i++)
      cout << std::setw(8) << std::setprecision(5) << vec_d[a][1][i] << "\t ";
    cout << endl;
    cout << "\t C : ";
    for (int i = 0; i <= Derivative_Max; i++)
      cout << std::setw(8) << std::setprecision(5) << vec_d[a][2][i] << "\t ";
    cout << endl;
    cout << "\t D : ";
    for (int i = 0; i <= Derivative_Max; i++)
      cout << std::setw(8) << std::setprecision(5) << vec_d[a][3][i] << "\t ";
    cout << endl;
    cout << "\t E : ";
    for (int i = 0; i <= Derivative_Max; i++)
      cout << std::setw(8) << std::setprecision(5) << vec_d[a][4][i] << "\t ";
    cout << endl;
  }
}

double MBPT_SEC::cal_2Bfr1B_normalized_Tz(int a, int b, int c, int d, int J2,
                                          vector<double> &V_2f1B) {

  int t2z_a = SP.state_Tz[a].tz_2;
  int t2z_b = SP.state_Tz[b].tz_2;

  int T2z = t2z_a + t2z_b;

  int j2_c = SP.state_Tz[c].j_2;
  int t2z_c = SP.state_Tz[c].tz_2;

  int j2_d = SP.state_Tz[d].j_2;
  int t2z_d = SP.state_Tz[d].tz_2;

  double res = 0.0;
  vector<double> V_dir(Derivative_Max + 1, 0.0);
  vector<double> V_ex(Derivative_Max + 1, 0.0);
  V_2f1B.clear();
  V_2f1B.resize(Derivative_Max + 1);

  double Norm_fac = 1;
  if (a == b)
    Norm_fac = Norm_fac / Norm_fac_C;
  if (c == d)
    Norm_fac = Norm_fac / Norm_fac_C;

  if (a == c) {
    V_dir = V_dir + Sm_vec[a];
  }
  if (b == d) {
    V_dir = V_dir + Sm_vec[b];
  }

  if (a == d) {
    V_ex = V_ex + Sm_vec[a];
  }

  if (b == c) {
    V_ex = V_ex + Sm_vec[b];
  }

  if (T2z != 0) {
    int j2c = SP.state_Tz[c].j_2;
    int j2d = SP.state_Tz[d].j_2;
    int phase_t = (j2c + j2d - J2) / 2;
    int phase_ex = phase_cal(phase_t);
    for (int i = 0; i <= Derivative_Max; i++) {
      V_2f1B[i] = Norm_fac * (V_dir[i] - phase_ex * V_ex[i]);
    }
  } else {
    for (int i = 0; i <= Derivative_Max; i++) {
      V_2f1B[i] = Norm_fac * (V_dir[i]);
    }
  }

  res = V_2f1B[0];

  return res;
}
void MBPT_SEC::build_sec_Qm() {
  if (Sm_build_flag != 1) {
    cout << "Wrong @ MBPT_SEC::build_sec_Qm : Build Sm first" << endl;
    exit(0);
  }
  int ch_size = DP_Model.channel_Tz.size();
  // cout << "ch_size = " << ch_size << endl;
  Qm_vec.resize(ch_size);
  Sm2Bf1B_vec.resize(ch_size);
  for (int ch = 0; ch < ch_size; ch++) {
    int T2z = DP_Model.channel_Tz[ch].T_2z;
    int J2 = DP_Model.channel_Tz[ch].J_2;
    int size = DP_Model.state_Tz[ch].size();
    // cout << "\t size = " << size << endl;
    Qm_vec[ch].resize(size);
    Sm2Bf1B_vec[ch].resize(size);
    for (int i = 0; i < size; i++) {
      Qm_vec[ch][i].resize(size);
      Sm2Bf1B_vec[ch][i].resize(size);
      //
      int a = DP_Model.state_Tz[ch][i].a;
      int b = DP_Model.state_Tz[ch][i].b;
      for (int j = 0; j < size; j++) {
        int c = DP_Model.state_Tz[ch][j].a;
        int d = DP_Model.state_Tz[ch][j].b;
        // double val_1st = get2B(a, b, c, d, J2, T2z);
        double val_1st = get2B_normalized_Tz(a, b, c, d, J2, T2z);
        vector<double> vec_d(Derivative_Max + 1, 0.0);
        vector<double> vec_d_2Bf1B(Derivative_Max + 1, 0.0);

        //
        Qm_vec[ch][i][j].resize(Derivative_Max + 1);
        Sm2Bf1B_vec[ch][i][j].resize(Derivative_Max + 1);
        //
        double val_Tw = cal_TwoBody_normalized_Tz(a, b, c, d, J2, 100, vec_d);
        double val_2Bf1B =
            cal_2Bfr1B_normalized_Tz(a, b, c, d, J2, vec_d_2Bf1B);
        vec_d[0] = val_1st + vec_d[0] + vec_d_2Bf1B[0];
        // vec_d[0] = vec_d[0];
        // Qm_vec[ch][i][j] = vec_d + vec_d_2Bf1B;
        for (int k = 0; k <= Derivative_Max; k++) {
          // Times m!
          Qm_vec[ch][i][j][k] = vec_d[k] * Frac_vec[k] + vec_d_2Bf1B[k];
          Sm2Bf1B_vec[ch][i][j][k] = vec_d_2Bf1B[k];
        }
      }
    }
  }
}
void MBPT_SEC::build_sec_TwB_sep() {
  int ch_size = DP_Model.channel_Tz.size();
  cout << "ch_size = " << ch_size << endl;
  vec_TwA_d.resize(ch_size);
  vec_TwB_d.resize(ch_size);
  vec_TwC_d.resize(ch_size);
  vec_TwD_d.resize(ch_size);
  vec_TwE_d.resize(ch_size);
  vec_TwF_d.resize(ch_size);
  for (int ch = 0; ch < ch_size; ch++) {
    int T2z = DP_Model.channel_Tz[ch].T_2z;
    int J2 = DP_Model.channel_Tz[ch].J_2;
    int size = DP_Model.state_Tz[ch].size();
    cout << "\t size = " << size << endl;
    vec_TwA_d[ch].resize(size);
    vec_TwB_d[ch].resize(size);
    vec_TwC_d[ch].resize(size);
    vec_TwD_d[ch].resize(size);
    vec_TwE_d[ch].resize(size);
    vec_TwF_d[ch].resize(size);
    for (int i = 0; i < size; i++) {

      vec_TwA_d[ch][i].resize(size);
      vec_TwB_d[ch][i].resize(size);
      vec_TwC_d[ch][i].resize(size);
      vec_TwD_d[ch][i].resize(size);
      vec_TwE_d[ch][i].resize(size);
      vec_TwF_d[ch][i].resize(size);
      //
      int a = DP_Model.state_Tz[ch][i].a;
      int b = DP_Model.state_Tz[ch][i].b;
      for (int j = 0; j < size; j++) {
        int c = DP_Model.state_Tz[ch][j].a;
        int d = DP_Model.state_Tz[ch][j].b;
        // double val_1st = get2B(a, b, c, d, J2, T2z);
        vector<vector<double>> vec_d(8);
        for (int k = 0; k < 8; k++) {
          vec_d[k].resize(Derivative_Max + 1);
          for (int l = 0; l <= Derivative_Max; l++)
            vec_d[k][l] = 0.0;
        }
        //
        vec_TwA_d[ch][i][j].resize(Derivative_Max + 1);
        vec_TwB_d[ch][i][j].resize(Derivative_Max + 1);
        vec_TwC_d[ch][i][j].resize(Derivative_Max + 1);
        vec_TwD_d[ch][i][j].resize(Derivative_Max + 1);
        vec_TwE_d[ch][i][j].resize(Derivative_Max + 1);
        vec_TwF_d[ch][i][j].resize(Derivative_Max + 1);
        //
        double val_TA = cal_TwoBody_normalized_Tz(a, b, c, d, J2, 0, vec_d[0]);
        double val_TB_1 =
            cal_TwoBody_normalized_Tz(a, b, c, d, J2, 1, vec_d[1]);
        double val_TB_2 =
            cal_TwoBody_normalized_Tz(a, b, c, d, J2, 2, vec_d[2]);
        double val_TC_1 =
            cal_TwoBody_normalized_Tz(a, b, c, d, J2, 3, vec_d[3]);
        double val_TC_2 =
            cal_TwoBody_normalized_Tz(a, b, c, d, J2, 4, vec_d[4]);
        double val_T_ph = val_TB_1 + val_TB_2 + val_TC_1 + val_TC_2;
        double val_TD = cal_TwoBody_normalized_Tz(a, b, c, d, J2, 5, vec_d[5]);
        double val_TE = cal_TwoBody_normalized_Tz(a, b, c, d, J2, 6, vec_d[6]);
        double val_TF = cal_TwoBody_normalized_Tz(a, b, c, d, J2, 7, vec_d[7]);
        vec_TwA_d[ch][i][j] = vec_d[0];
        vec_TwB_d[ch][i][j] = vec_d[1] + vec_d[2];
        vec_TwC_d[ch][i][j] = vec_d[3] + vec_d[4];
        vec_TwD_d[ch][i][j] = vec_d[5];
        vec_TwE_d[ch][i][j] = vec_d[6];
        vec_TwF_d[ch][i][j] = vec_d[7];
      }
    }
  }
}
double MBPT_SEC::cal_OneBody(int a, vector<double> &vec_d) {
  int sp_size = SP.state_Tz.size();
  double res = 0.0;
  vec_d.clear();
  vec_d.resize(Derivative_Max + 1);
  vector<vector<double>> vec_sep;

  vec_sep.resize(5);
  for (int k = 0; k < 5; k++) {
    vec_sep[k].resize(Derivative_Max + 1);
    for (int j = 0; j <= Derivative_Max; j++)
      vec_sep[k][j] = 0.0;
  }

  vector<double> vec_res(5, 0.0);
  double E_start = SP.state_Tz[a].E * hbar_omega;
  vec_res[0] = cal_OA(a, E_start, vec_sep[0]);
  vec_res[1] = cal_OB(a, E_start, vec_sep[1]);
  vec_res[2] = cal_OC(a, E_start, vec_sep[2]);
  vec_res[3] = cal_OD(a, E_start, vec_sep[3]);
  vec_res[4] = cal_OE(a, E_start, vec_sep[4]);

  for (int j = 0; j <= Derivative_Max; j++) {
    res = vec_res[0] + vec_res[1] + vec_res[2] + vec_res[3] + vec_res[4];
    vec_d[j] = vec_sep[0][j] + vec_sep[1][j] + vec_sep[2][j] + vec_sep[3][j] +
               vec_sep[4][j];
  }
  return res;
}
void MBPT_SEC::print_sec_OB_deriva() {
  cout << " this->hbar_omega = " << this->hbar_omega << endl;
  int sp_size = SP.state_Tz.size();

  // cout << "a\t A\t\t B\t\t C\t\t D\t\t  E*2\t " << endl;
  for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
    cout << " --- SP Tz : " << a << "\t ---" << endl;
    cout << "\t S-box : ";
    for (int i = 0; i <= Derivative_Max; i++)
      cout << std::setw(8) << std::setprecision(5) << Sm_vec[a][i] << "\t ";
    cout << endl;
  }
}
double MBPT_SEC::cal_OA(int a, double E_start, vector<double> &vec_d) {
  int phase = 1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;
  // double res_t = 0.0;
  for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
    // cout<<"\t -- p : "<<p<<"\t -- "<<endl;
    int e_p = SP.state_Tz[p].E;
    int j2_p = SP.state_Tz[p].j_2;
    int t2z_p = SP.state_Tz[p].tz_2;
    double E_mid = this->hbar_omega * (e_p);
    double E_dom = E_start - E_mid;
    if (E_dom == 0.0)
      continue;
    if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
      continue;
    if (E_dom == 0.0)
      continue;
    if (j2_a != j2_p || t2z_a != t2z_p)
      continue;
    // int flag = 0;
    double val_U = get1B(a, p);
    double res_add = val_U * val_U;
    res_add = res_add / E_dom;
    res += res_add;
    //
    vec_d[0] += res_add;
    int Derive_phase = 1;
    double Derive_Edom = 1.0;
    // cout << " +++ E_dom : " << E_dom << "Derive_Edom = " << Derive_Edom
    //      << "\t ";
    for (int i = 1; i <= Derivative_Max; i++) {
      Derive_phase = Derive_phase * -1;
      Derive_Edom = Derive_Edom / E_dom;
      // cout << Derive_Edom << "\t ";
      vec_d[i] += res_add * Derive_phase * Derive_Edom;
    }
    // cout << endl;
  }
  double result = res * phase;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = phase * vec_d[i];
  }
  return result;
}
double MBPT_SEC::cal_OB(int a, double E_start, vector<double> &vec_d) {
  int phase = -1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;
  // double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    double E_mid = this->hbar_omega * (e_a + e_a - e_h);
    double E_dom = E_start - E_mid;
    if (E_dom == 0.0)
      continue;
    if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
      continue;
    if (j2_a != j2_h || t2z_a != t2z_h)
      continue;
    double val_U = get1B(a, h);
    double res_add = val_U * val_U;
    res_add = res_add / E_dom;
    res += res_add;
    //
    vec_d[0] += res_add;
    int Derive_phase = 1;
    double Derive_Edom = 1.0;
    for (int i = 1; i <= Derivative_Max; i++) {
      Derive_phase = Derive_phase * -1;
      Derive_Edom = Derive_Edom / E_dom;
      vec_d[i] += res_add * Derive_phase * Derive_Edom;
    }
  }
  double result = res * phase;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = phase * vec_d[i];
  }
  return result;
}

double MBPT_SEC::cal_OC(int a, double E_start, vector<double> &vec_d) {
  int phase = 1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  double res = 0.0;
  // double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int J2_min = std::abs(j2_a - j2_h);
    int J2_max = j2_a + j2_h;
    int T2z = t2z_a + t2z_h;
    int par_t = SP.state_Tz[a].l + SP.state_Tz[h].l;
    int Par_x = par_t % 2;

    for (int p1 = sp_P_beg_Tz; p1 <= sp_Q_end_Tz; p1++) {
      int e_p1 = SP.state_Tz[p1].E;
      int j2_p1 = SP.state_Tz[p1].j_2;
      int t2z_p1 = SP.state_Tz[p1].tz_2;
      for (int p2 = sp_P_beg_Tz; p2 <= sp_Q_end_Tz; p2++) {
        int e_p2 = SP.state_Tz[p2].E;
        int j2_p2 = SP.state_Tz[p2].j_2;
        int t2z_p2 = SP.state_Tz[p2].tz_2;
        if (T2z != (t2z_p1 + t2z_p2))
          continue;

        par_t = SP.state_Tz[p1].l + SP.state_Tz[p2].l;
        int Par_y = par_t % 2;
        if (Par_x != Par_y)
          continue;
        double E_mid = this->hbar_omega * (e_p1 + e_p2 - e_h);
        double E_dom = E_start - E_mid;
        if (E_dom == 0.0)
          continue;
        if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
          continue;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          // cout << "-- C1 --  J2 : " << J2 << "\t j2_p1 : " << j2_p1
          //      << "\t j2_p2 : " << j2_p2 << "\t a : " << a << "\t h : " <<
          //      h
          //      << "\t p1 : " << p1 << "\t p2 : " << p2 << endl;
          if (Tri(J2, j2_p1, j2_p2) < 0)
            continue;

          double val_2V = get2B(a, h, p1, p2, J2, T2z);

          double res_t = val_2V * val_2V;
          double res_add = (J2 + 1) * res_t / E_dom;
          res += res_add;
          // cout << "\t C1 += " << res_t << "\t";
          // cout << "a : " << a << "\t h : " << h << "\t p1 : " << p1
          //      << "\t p2 : " << p2 << endl;
          //
          vec_d[0] += res_add;
          int Derive_phase = 1;
          double Derive_Edom = 1.0;
          for (int i = 1; i <= Derivative_Max; i++) {
            Derive_phase = Derive_phase * -1;
            Derive_Edom = Derive_Edom / E_dom;
            vec_d[i] += res_add * Derive_phase * Derive_Edom;
          }
        }
      }
    }
  }
  double pre_fac = 1.0 / (2.0 * (j2_a + 1.0));
  double result = res * phase * pre_fac;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = phase * pre_fac * vec_d[i];
  }
  return result;
}

double MBPT_SEC::cal_OD(int a, double E_start, vector<double> &vec_d) {
  int phase = -1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  double res = 0.0;
  // double res_t = 0.0;
  for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
    // cout << " -- p:" << p << endl;
    int e_p = SP.state_Tz[p].E;
    int j2_p = SP.state_Tz[p].j_2;
    int t2z_p = SP.state_Tz[p].tz_2;
    int J2_min = std::abs(j2_a - j2_p);
    int J2_max = j2_a + j2_p;
    int T2z = t2z_a + t2z_p;
    int par_t = SP.state_Tz[a].l + SP.state_Tz[p].l;
    int Par_x = par_t % 2;

    for (int h1 = 0; h1 <= sp_core_Tz; h1++) {
      int e_h1 = SP.state_Tz[h1].E;
      int j2_h1 = SP.state_Tz[h1].j_2;
      int t2z_h1 = SP.state_Tz[h1].tz_2;
      for (int h2 = 0; h2 <= sp_core_Tz; h2++) {
        int e_h2 = SP.state_Tz[h2].E;
        int j2_h2 = SP.state_Tz[h2].j_2;
        int t2z_h2 = SP.state_Tz[h2].tz_2;
        double E_mid = this->hbar_omega * (e_a + e_a - e_h1 - e_h2 + e_p);
        double E_dom = E_start - E_mid;
        if (E_dom == 0.0)
          continue;
        if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
          continue;
        if (T2z != (t2z_h1 + t2z_h2))
          continue;
        par_t = SP.state_Tz[h1].l + SP.state_Tz[h2].l;
        int Par_y = par_t % 2;
        if (Par_x != Par_y)
          continue;
        // cout << " -- h1:" << h1 << " h2:" << h2 << endl;
        for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
          // cout << "-- C2 --  J2 : " << J2 << "\t j2_h1 : " << j2_h1
          //      << "\t j2_h1 : " << j2_h2 << "\t a : " << a << "\t p : " <<
          //      p
          //      << "\t h1 : " << h1 << "\t h2 : " << h2 << endl;
          if (Tri(J2, j2_h1, j2_h2) < 0)
            continue;

          double val_2V = get2B(h1, h2, a, p, J2, T2z);

          double res_t = val_2V * val_2V;
          double res_add = (J2 + 1) * res_t / E_dom;
          res += res_add;
          // cout << " ------ ------- " << endl;
          // cout << "\t T2z = " << T2z << "  J2 = " << J2 << " a : " << a
          //      << "\t p : " << p << "\t h1 : " << h1 << "\t h2 : " << h2
          //      << endl;
          // cout << "\t C2 += " << res_t << "\t"
          //      << "vec_2V2N = " << vec_2V2N[ch_DP][index_y][index_x]
          //      << "   vec_2V3N = " << vec_2V3N[ch_DP][index_y][index_x] <<
          //      endl;
          //
          vec_d[0] += res_add;
          int Derive_phase = 1;
          double Derive_Edom = 1.0;
          for (int i = 1; i <= Derivative_Max; i++) {
            Derive_phase = Derive_phase * -1;
            Derive_Edom = Derive_Edom / E_dom;
            vec_d[i] += res_add * Derive_phase * Derive_Edom;
          }
        }
      }
    }
  }
  double pre_fac = 1.0 / (2.0 * (j2_a + 1.0));
  double result = res * phase * pre_fac;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = phase * pre_fac * vec_d[i];
  }
  return result;
}

double MBPT_SEC::cal_OE(int a, double E_start, vector<double> &vec_d) {
  int phase = 1;
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  double res = 0.0;
  // double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int T2z_ah = t2z_a + t2z_h;
    int J2_max = j2_a + j2_h;
    int J2_min = std::abs(j2_a - j2_h);
    int Par_h = SP.state_Tz[h].l % 2;
    for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
      int e_p = SP.state_Tz[p].E;
      int j2_p = SP.state_Tz[p].j_2;
      int t2z_p = SP.state_Tz[p].tz_2;
      int Par_p = SP.state_Tz[p].l % 2;
      double E_mid = this->hbar_omega * (e_a + e_p - e_h);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > (this->hbar_omega * this->E_cut))
        continue;
      int T2z_ap = t2z_a + t2z_p;
      if (j2_h != j2_p || t2z_h != t2z_p || T2z_ap != T2z_ah)
        continue;
      if (Par_h != Par_p)
        continue;

      double val_U = get1B(h, p);

      // double val_t = 0.0;
      for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
        if (Tri(J2, j2_a, j2_p) < 0)
          continue;

        double val_2V = get2B(a, h, a, p, J2, T2z_ap);

        double res_t = val_2V * val_U;
        double res_add = (J2 + 1) * res_t / E_dom;
        //
        // cout << "------  ------" << endl;
        // cout << "  J2=" << J2 << "\t a_z = " << a << "  p_z = " << p
        //      << "  h_z = " << h << endl;
        // cout << "\t val_U = " << val_U
        //      << "\t vec_2V3N = " << vec_2V3N[ch_DP][index_y][index_x] <<
        //      endl;
        res += res_add;
        //
        vec_d[0] += res_add;
        int Derive_phase = 1;
        double Derive_Edom = 1.0;
        for (int i = 1; i <= Derivative_Max; i++) {
          Derive_phase = Derive_phase * -1;
          Derive_Edom = Derive_Edom / E_dom;
          vec_d[i] += res_add * Derive_phase * Derive_Edom;
        }
      }
    }
  }
  // the exchange of two vertices !!! *2
  double pre_fac = 2 * 1.0 / (j2_a + 1.0);
  double result = res * phase * pre_fac;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = phase * pre_fac * vec_d[i];
  }
  return result;
}

double MBPT_SEC::cal_TwoBody(int a, int b, int c, int d, int J2, int Flag,
                             vector<double> &vec_d) {
  double res = 0.0;
  vector<vector<double>> vec_sep;
  vec_sep.resize(8);
  for (int i = 0; i < 8; i++) {
    vec_sep[i].resize(Derivative_Max + 1);
    for (int j = 0; j <= Derivative_Max; j++)
      vec_sep[i][j] = 0.0;
  }
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;
  int e_b = SP.state_Tz[b].E;
  int j2_b = SP.state_Tz[b].j_2;
  int t2z_b = SP.state_Tz[b].tz_2;
  int T2z = t2z_a + t2z_b;
  int Par = (SP.state_Tz[a].l + SP.state_Tz[b].l) % 2;
  int e_c = SP.state_Tz[c].E;
  int j2_c = SP.state_Tz[c].j_2;
  int t2z_c = SP.state_Tz[c].tz_2;
  int e_d = SP.state_Tz[d].E;
  int j2_d = SP.state_Tz[d].j_2;
  int t2z_d = SP.state_Tz[d].tz_2;

  double E_start = (e_c + e_d) * this->hbar_omega;

  if (Par != (SP.state_Tz[c].l + SP.state_Tz[d].l) % 2)
    return res;

  if (T2z != (t2z_c + t2z_d))
    return res;

  if (Tri(j2_a, j2_b, J2) < 0)
    return res;
  if (Tri(j2_c, j2_d, J2) < 0)
    return res;
  vector<double> res_vec(9, 0);
  switch (Flag) {
  case 0:
    res_vec[0] = cal_TwA(a, b, c, d, J2, T2z, Par, E_start, vec_d);
    break;
  case 1:
    res_vec[1] = cal_TwB_1(a, b, c, d, J2, T2z, Par, E_start, vec_d);
    break;
  case 2:
    res_vec[2] = cal_TwB_2(a, b, c, d, J2, T2z, Par, E_start, vec_d);
    break;
  case 3:
    res_vec[3] = cal_TwC_1(a, b, c, d, J2, T2z, Par, E_start, vec_d);
    break;
  case 4:
    res_vec[4] = cal_TwC_2(a, b, c, d, J2, T2z, Par, E_start, vec_d);
    break;
  case 5:
    res_vec[5] = cal_TwD(a, b, c, d, J2, T2z, Par, E_start, vec_d);
    break;
  case 6:
    res_vec[6] = cal_TwE(a, b, c, d, J2, T2z, Par, E_start, vec_d);
    break;
  case 7:
    res_vec[7] = cal_TwF(a, b, c, d, J2, T2z, Par, E_start, vec_d);
    break;
  case 100:

    res_vec[0] = cal_TwA(a, b, c, d, J2, T2z, Par, E_start, vec_sep[0]);
    res_vec[1] = cal_TwB_1(a, b, c, d, J2, T2z, Par, E_start, vec_sep[1]);
    res_vec[2] = cal_TwB_2(a, b, c, d, J2, T2z, Par, E_start, vec_sep[2]);
    res_vec[3] = cal_TwC_1(a, b, c, d, J2, T2z, Par, E_start, vec_sep[3]);
    res_vec[4] = cal_TwC_2(a, b, c, d, J2, T2z, Par, E_start, vec_sep[4]);
    res_vec[5] = cal_TwD(a, b, c, d, J2, T2z, Par, E_start, vec_sep[5]);
    res_vec[6] = cal_TwE(a, b, c, d, J2, T2z, Par, E_start, vec_sep[6]);
    res_vec[7] = cal_TwF(a, b, c, d, J2, T2z, Par, E_start, vec_sep[7]);

    vec_d.clear();
    vec_d.resize(Derivative_Max + 1);
    for (int j = 0; j <= Derivative_Max; j++) {
      for (int i = 0; i < 8; i++) {
        vec_d[j] += vec_sep[i][j];
      }
    }

    break;
  default:
    cout << "Wrong @ MBPT_SEC::cal_TwoBody " << endl;
    cout << "Flag = " << Flag << endl;
    exit(0);
    break;
  }
  for (int i = 0; i < 8; i++)
    res += res_vec[i];
  return res;
}

double MBPT_SEC::cal_TwoBody_normalized_Tz(int a, int b, int c, int d, int J2,
                                           int Flag, vector<double> &vec_d) {
  double tbme = this->cal_TwoBody(a, b, c, d, J2, Flag, vec_d);
  double Norm_fac = 1.0;
  double Norm_fac_t = 1.0;
  if (a == b)
    Norm_fac_t = Norm_fac_t * (std::sqrt(2));
  if (c == d)
    Norm_fac_t = Norm_fac_t * (std::sqrt(2));
  Norm_fac = 1.0 / (Norm_fac_t);
  double tbme_norm = tbme * Norm_fac;
  for (int i = 0; i <= Derivative_Max; i++) {
    if (std::abs(vec_d[i]) < EPS)
      vec_d[i] = 0.0;
    vec_d[i] = vec_d[i] * Norm_fac;
  }
  return tbme_norm;
}

double MBPT_SEC::cal_TwA(int a, int b, int c, int d, int J2, int T2z, int Par,
                         double E_start, vector<double> &vec_d) {
  double res = 0.0;
  // cout<<"Start cal_TwA"<<endl;
  // cout << "vec_d.size() = " << vec_d.size() << endl;
  // vec_d.clear();
  // vec_d.resize(Derivative_Max + 1);
  for (int p1 = sp_P_beg_Tz; p1 <= sp_Q_end_Tz; p1++) {
    int e_p1 = SP.state_Tz[p1].E;
    int j2_p1 = SP.state_Tz[p1].j_2;
    int t2z_p1 = SP.state_Tz[p1].tz_2;
    for (int p2 = sp_P_beg_Tz; p2 <= sp_Q_end_Tz; p2++) {
      int e_p2 = SP.state_Tz[p2].E;
      int j2_p2 = SP.state_Tz[p2].j_2;
      int t2z_p2 = SP.state_Tz[p2].tz_2;
      int Par_p1p2 = (SP.state_Tz[p1].l + SP.state_Tz[p2].l) % 2;
      if (Par != Par_p1p2)
        continue;
      double E_mid = this->hbar_omega * (e_p1 + e_p2);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      double mat1 = get2B(a, b, p1, p2, J2, T2z);
      double mat2 = get2B(c, d, p1, p2, J2, T2z);
      double res_add = mat1 * mat2 / E_dom;
      res += res_add;
      vec_d[0] += res_add;
      int Derive_phase = 1;
      double Derive_Edom = 1.0;
      for (int i = 1; i <= Derivative_Max; i++) {
        Derive_phase = Derive_phase * -1;
        Derive_Edom = Derive_Edom / E_dom;
        vec_d[i] += res_add * Derive_phase * Derive_Edom;
      }
    }
  }
  res = 0.5 * res;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = 0.5 * vec_d[i];
  }
  // cout << "calculate TA res = " << res << endl;
  return res;
}
double MBPT_SEC::cal_TwB_1(int a, int b, int c, int d, int J2, int T2z, int Par,
                           double E_start, vector<double> &vec_d) {

  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  int e_b = SP.state_Tz[b].E;
  int j2_b = SP.state_Tz[b].j_2;
  int t2z_b = SP.state_Tz[b].tz_2;

  int e_c = SP.state_Tz[c].E;
  int j2_c = SP.state_Tz[c].j_2;
  int t2z_c = SP.state_Tz[c].tz_2;

  int e_d = SP.state_Tz[d].E;
  int j2_d = SP.state_Tz[d].j_2;
  int t2z_d = SP.state_Tz[d].tz_2;
  //

  double res = 0.0;
  int phase_1 = -1;
  double pre_fac = 1.0 / Sqrt_vec[J2 + 1];
  //
  int J2min_ac = std::abs(j2_a - j2_c);
  int J2min_bd = std::abs(j2_b - j2_d);
  int J2min = std::max(J2min_ac, J2min_bd);

  int J2max_ac = (j2_a + j2_c);
  int J2max_bd = (j2_b + j2_d);
  int J2max = std::min(J2max_ac, J2max_bd);
  // cout << " +++ J2min = " << J2min << "\t J2max = " << J2max << endl;
  if (J2min > J2max) {
    cout << "Wrong @  MBPT_SEC::cal_TB_1" << endl;
    exit(0);
  }
  //
  double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int T2z_ah = t2z_a + t2z_h;
    int T2z_hd = t2z_d + t2z_h;

    for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
      int e_p = SP.state_Tz[p].E;
      int j2_p = SP.state_Tz[p].j_2;
      int t2z_p = SP.state_Tz[p].tz_2;
      if (T2z_ah != t2z_c + t2z_p)
        continue;
      if (T2z_hd != t2z_b + t2z_p)
        continue;
      double E_mid = this->hbar_omega * (e_c + e_p - e_h + e_b);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      for (int J2p = J2min; J2p <= J2max; J2p = J2p + 2) {
        if (Tri(J2p, j2_h, j2_p) < 0)
          continue;
        double val_X0 = cal_X0(j2_c, j2_a, J2p, j2_d, j2_b, J2p, J2, J2);
        // double val_X0 = cal_X0(j2_c, j2_d, J2, j2_a, j2_b, J2, J2p, J2p);
        if (val_X0 == 0.0)
          continue;
        double mat1 = cal_CRC_ME_I(a, h, c, p, J2p, T2z_ah);
        double mat2 = cal_CRC_ME_I(p, b, h, d, J2p, T2z_hd);
        //
        int phase_t = (j2_p + j2_h - J2p) / 2;
        int phase_2 = phase_cal(phase_t);
        double res_add =
            phase_2 * Sqrt_vec[J2p + 1] * val_X0 * mat1 * mat2 / E_dom;
        res_t += res_add;
        // cout << " ---- J2:" << J2 << "\t a:" << a << " b:" << b << " c:" << c
        //      << " d:" << d << " ---- " << endl;
        // cout << "\t h:" << h << " p:" << p << "\t J2p:" << J2p << endl;
        // cout << "\t phase_2:" << phase_2
        //      << "\t Sqrt_vec[J2p+1]:" << Sqrt_vec[J2p + 1]
        //      << "\t E_dom = " << E_dom << "\t val_X0:" << val_X0
        //      << "\t mat1:" << mat1 << "\t mat2:" << mat2
        //      << "\t res_add = " << res_add << endl;
        vec_d[0] += res_add;
        int Derive_phase = 1;
        double Derive_Edom = 1.0;
        for (int i = 1; i <= Derivative_Max; i++) {
          Derive_phase = Derive_phase * -1;
          Derive_Edom = Derive_Edom / E_dom;
          vec_d[i] += res_add * Derive_phase * Derive_Edom;
        }
      }
    }
  }
  res = phase_1 * pre_fac * res_t;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = phase_1 * pre_fac * vec_d[i];
  }
  // cout << "TB_1 = " << res << endl;
  return res;
}
double MBPT_SEC::cal_TwB_2(int a, int b, int c, int d, int J2, int T2z, int Par,
                           double E_start, vector<double> &vec_d) {
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  int e_b = SP.state_Tz[b].E;
  int j2_b = SP.state_Tz[b].j_2;
  int t2z_b = SP.state_Tz[b].tz_2;

  int e_c = SP.state_Tz[c].E;
  int j2_c = SP.state_Tz[c].j_2;
  int t2z_c = SP.state_Tz[c].tz_2;

  int e_d = SP.state_Tz[d].E;
  int j2_d = SP.state_Tz[d].j_2;
  int t2z_d = SP.state_Tz[d].tz_2;
  //

  double res = 0.0;
  int phase_1t = (j2_c + j2_d - J2) / 2;
  int phase_1 = phase_cal(phase_1t);
  double pre_fac = 1.0 / Sqrt_vec[J2 + 1];
  //
  int J2min_ad = std::abs(j2_a - j2_d);
  int J2min_bc = std::abs(j2_b - j2_c);
  int J2min = std::max(J2min_ad, J2min_bc);

  int J2max_ad = (j2_a + j2_d);
  int J2max_bc = (j2_b + j2_c);
  int J2max = std::min(J2max_ad, J2max_bc);
  if (J2min > J2max) {
    cout << "Wrong @  MBPT_SEC::cal_TB_2" << endl;
    exit(0);
  }
  //
  double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int T2z_ah = t2z_a + t2z_h;
    int T2z_hc = t2z_c + t2z_h;

    for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
      int e_p = SP.state_Tz[p].E;
      int j2_p = SP.state_Tz[p].j_2;
      int t2z_p = SP.state_Tz[p].tz_2;
      if (T2z_ah != t2z_d + t2z_p)
        continue;
      if (T2z_hc != t2z_b + t2z_p)
        continue;
      double E_mid = this->hbar_omega * (e_d + e_p - e_h + e_b);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      for (int J2p = J2min; J2p <= J2max; J2p = J2p + 2) {
        if (Tri(J2p, j2_h, j2_p) < 0)
          continue;
        double val_X0 = cal_X0(j2_d, j2_a, J2p, j2_c, j2_b, J2p, J2, J2);
        if (val_X0 == 0.0)
          continue;
        double mat1 = cal_CRC_ME_I(a, h, d, p, J2p, T2z_ah);
        double mat2 = cal_CRC_ME_I(p, b, h, c, J2p, T2z_hc);
        //
        int phase_t = (j2_p + j2_h - J2p) / 2;
        int phase_2 = phase_cal(phase_t);
        double res_add =
            phase_2 * Sqrt_vec[J2p + 1] * val_X0 * mat1 * mat2 / E_dom;
        res_t += res_add;
        //
        vec_d[0] += res_add;
        int Derive_phase = 1;
        double Derive_Edom = 1.0;
        for (int i = 1; i <= Derivative_Max; i++) {
          Derive_phase = Derive_phase * -1;
          Derive_Edom = Derive_Edom / E_dom;
          vec_d[i] += res_add * Derive_phase * Derive_Edom;
        }
      }
    }
  }
  res = phase_1 * pre_fac * res_t;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = phase_1 * pre_fac * vec_d[i];
  }
  // cout << "TB_2 = " << res << endl;
  return res;
}
double MBPT_SEC::cal_TwC_1(int a, int b, int c, int d, int J2, int T2z, int Par,
                           double E_start, vector<double> &vec_d) {
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  int e_b = SP.state_Tz[b].E;
  int j2_b = SP.state_Tz[b].j_2;
  int t2z_b = SP.state_Tz[b].tz_2;

  int e_c = SP.state_Tz[c].E;
  int j2_c = SP.state_Tz[c].j_2;
  int t2z_c = SP.state_Tz[c].tz_2;

  int e_d = SP.state_Tz[d].E;
  int j2_d = SP.state_Tz[d].j_2;
  int t2z_d = SP.state_Tz[d].tz_2;
  //

  double res = 0.0;
  int phase_1 = -1;
  double pre_fac = 1.0 / Sqrt_vec[J2 + 1];
  //
  int J2min_ac = std::abs(j2_a - j2_c);
  int J2min_bd = std::abs(j2_b - j2_d);
  int J2min = std::max(J2min_ac, J2min_bd);

  int J2max_ac = std::abs(j2_a + j2_c);
  int J2max_bd = std::abs(j2_b + j2_d);
  int J2max = std::min(J2max_ac, J2max_bd);
  if (J2min > J2max) {
    cout << "Wrong @  MBPT_SEC::cal_TC_1" << endl;
    exit(0);
  }
  //
  double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int T2z_ch = t2z_c + t2z_h;
    int T2z_hb = t2z_b + t2z_h;

    for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
      int e_p = SP.state_Tz[p].E;
      int j2_p = SP.state_Tz[p].j_2;
      int t2z_p = SP.state_Tz[p].tz_2;
      if (T2z_ch != t2z_a + t2z_p)
        continue;
      if (T2z_hb != t2z_d + t2z_p)
        continue;
      double E_mid = this->hbar_omega * (e_a + e_p - e_h + e_d);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      for (int J2p = J2min; J2p <= J2max; J2p = J2p + 2) {
        if (Tri(J2p, j2_h, j2_p) < 0)
          continue;
        double val_X0 = cal_X0(j2_c, j2_a, J2p, j2_d, j2_b, J2p, J2, J2);
        if (val_X0 == 0.0)
          continue;
        double mat1 = cal_CRC_ME_I(a, p, c, h, J2p, T2z_ch);
        double mat2 = cal_CRC_ME_I(h, b, p, d, J2p, T2z_hb);
        //
        int phase_t = (j2_p + j2_h - J2p) / 2;
        int phase_2 = phase_cal(phase_t);
        double res_add =
            phase_2 * Sqrt_vec[J2p + 1] * val_X0 * mat1 * mat2 / E_dom;
        res_t += res_add;
        // cout << " ---- J2:" << J2 << "\t a:" << a << " b:" << b << " c:" << c
        //      << " d:" << d << " ---- " << endl;
        // cout << "\t h:" << h << " p:" << p << "\t J2p:" << J2p << endl;
        // cout << "\t phase_2:" << phase_2
        //      << "\t Sqrt_vec[J2p+1]:" << Sqrt_vec[J2p + 1]
        //      << "\t E_dom = " << E_dom << "\t val_X0:" << val_X0
        //      << "\t mat1:" << mat1 << "\t mat2:" << mat2
        //      << "\t res_add = " << res_add << endl;
        //
        vec_d[0] += res_add;
        int Derive_phase = 1;
        double Derive_Edom = 1.0;
        for (int i = 1; i <= Derivative_Max; i++) {
          Derive_phase = Derive_phase * -1;
          Derive_Edom = Derive_Edom / E_dom;
          vec_d[i] += res_add * Derive_phase * Derive_Edom;
        }
      }
    }
  }
  res = phase_1 * pre_fac * res_t;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = phase_1 * pre_fac * vec_d[i];
  }
  // cout << "TC_1 = " << res << endl;
  return res;
}
double MBPT_SEC::cal_TwC_2(int a, int b, int c, int d, int J2, int T2z, int Par,
                           double E_start, vector<double> &vec_d) {
  int e_a = SP.state_Tz[a].E;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  int e_b = SP.state_Tz[b].E;
  int j2_b = SP.state_Tz[b].j_2;
  int t2z_b = SP.state_Tz[b].tz_2;

  int e_c = SP.state_Tz[c].E;
  int j2_c = SP.state_Tz[c].j_2;
  int t2z_c = SP.state_Tz[c].tz_2;

  int e_d = SP.state_Tz[d].E;
  int j2_d = SP.state_Tz[d].j_2;
  int t2z_d = SP.state_Tz[d].tz_2;
  //

  double res = 0.0;
  int phase_1t = (j2_c + j2_d - J2) / 2;
  int phase_1 = phase_cal(phase_1t);
  double pre_fac = 1.0 / Sqrt_vec[J2 + 1];
  //
  int J2min_ad = std::abs(j2_a - j2_d);
  int J2min_bc = std::abs(j2_b - j2_c);
  int J2min = std::max(J2min_ad, J2min_bc);

  int J2max_ad = std::abs(j2_a + j2_d);
  int J2max_bc = std::abs(j2_b + j2_c);
  int J2max = std::min(J2max_ad, J2max_bc);
  if (J2min > J2max) {
    cout << "Wrong @  MBPT_SEC::cal_TC_2" << endl;
    exit(0);
  }
  //
  double res_t = 0.0;
  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;
    int T2z_dh = t2z_d + t2z_h;
    int T2z_hb = t2z_b + t2z_h;

    for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
      int e_p = SP.state_Tz[p].E;
      int j2_p = SP.state_Tz[p].j_2;
      int t2z_p = SP.state_Tz[p].tz_2;
      if (T2z_dh != t2z_a + t2z_p)
        continue;
      if (T2z_hb != t2z_c + t2z_p)
        continue;
      double E_mid = this->hbar_omega * (e_a + e_p - e_h + e_c);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      for (int J2p = J2min; J2p <= J2max; J2p = J2p + 2) {
        if (Tri(J2p, j2_h, j2_p) < 0)
          continue;
        double val_X0 = cal_X0(j2_d, j2_a, J2p, j2_c, j2_b, J2p, J2, J2);
        if (val_X0 == 0.0)
          continue;
        double mat1 = cal_CRC_ME_I(a, p, d, h, J2p, T2z_dh);
        double mat2 = cal_CRC_ME_I(h, b, p, c, J2p, T2z_hb);
        //
        int phase_t = (j2_p + j2_h - J2p) / 2;
        int phase_2 = phase_cal(phase_t);
        double res_add =
            phase_2 * Sqrt_vec[J2p + 1] * val_X0 * mat1 * mat2 / E_dom;
        res_t += res_add;
        //
        vec_d[0] += res_add;
        int Derive_phase = 1;
        double Derive_Edom = 1.0;
        for (int i = 1; i <= Derivative_Max; i++) {
          Derive_phase = Derive_phase * -1;
          Derive_Edom = Derive_Edom / E_dom;
          vec_d[i] += res_add * Derive_phase * Derive_Edom;
        }
      }
    }
  }
  res = phase_1 * pre_fac * res_t;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = phase_1 * pre_fac * vec_d[i];
  }
  // cout << "TC_2 = " << res << endl;
  return res;
}
double MBPT_SEC::cal_TwD(int a, int b, int c, int d, int J2, int T2z, int Par,
                         double E_start, vector<double> &vec_d) {
  double res = 0.0;
  int e_a = SP.state_Tz[a].E;
  int e_b = SP.state_Tz[b].E;
  int e_c = SP.state_Tz[c].E;
  int e_d = SP.state_Tz[d].E;
  for (int h1 = 0; h1 <= sp_core_Tz; h1++) {
    int e_h1 = SP.state_Tz[h1].E;
    int j2_h1 = SP.state_Tz[h1].j_2;
    int t2z_h1 = SP.state_Tz[h1].tz_2;
    for (int h2 = 0; h2 <= sp_core_Tz; h2++) {
      int e_h2 = SP.state_Tz[h2].E;
      int j2_h2 = SP.state_Tz[h2].j_2;
      int t2z_h2 = SP.state_Tz[h2].tz_2;
      int Par_h1h2 = (SP.state_Tz[h1].l + SP.state_Tz[h2].l) % 2;
      if (Par != Par_h1h2)
        continue;
      double E_mid = this->hbar_omega * (e_a + e_b + e_c + e_d - e_h1 - e_h2);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      double mat1 = get2B(h1, h2, a, b, J2, T2z);
      double mat2 = get2B(h1, h2, c, d, J2, T2z);
      double res_add = mat1 * mat2 / E_dom;
      res += res_add;
      //
      vec_d[0] += res_add;
      int Derive_phase = 1;
      double Derive_Edom = 1.0;
      for (int i = 1; i <= Derivative_Max; i++) {
        Derive_phase = Derive_phase * -1;
        Derive_Edom = Derive_Edom / E_dom;
        vec_d[i] += res_add * Derive_phase * Derive_Edom;
      }
    }
  }
  res = 0.5 * res;
  for (int i = 0; i <= Derivative_Max; i++) {
    vec_d[i] = 0.5 * vec_d[i];
  }
  return res;
}
double MBPT_SEC::cal_TwE(int a, int b, int c, int d, int J2, int T2z, int Par,
                         double E_start, vector<double> &vec_d) {
  double res = 0.0;
  int e_a = SP.state_Tz[a].E;
  int l_a = SP.state_Tz[a].l;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  int e_b = SP.state_Tz[b].E;
  int l_b = SP.state_Tz[b].l;
  int j2_b = SP.state_Tz[b].j_2;
  int t2z_b = SP.state_Tz[b].tz_2;

  int e_c = SP.state_Tz[c].E;
  int l_c = SP.state_Tz[c].l;
  int j2_c = SP.state_Tz[c].j_2;
  int t2z_c = SP.state_Tz[c].tz_2;

  int e_d = SP.state_Tz[d].E;
  int l_d = SP.state_Tz[d].l;
  int j2_d = SP.state_Tz[d].j_2;
  int t2z_d = SP.state_Tz[d].tz_2;

  double res_E1 = 0.0;
  double res_E2 = 0.0;
  double res_E3 = 0.0;
  double res_E4 = 0.0;
  // vector<vector<double>> vec_sep(4);

  for (int p = sp_P_beg_Tz; p <= sp_Q_end_Tz; p++) {
    int e_p = SP.state_Tz[p].E;
    int l_p = SP.state_Tz[p].l;
    int j2_p = SP.state_Tz[p].j_2;
    int t2z_p = SP.state_Tz[p].tz_2;

    // E1
    if ((t2z_b == t2z_p) && (j2_b == j2_p) && (l_b % 2) == (l_p % 2)) {
      double E_mid = this->hbar_omega * (e_a + e_p);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      double mat1 = get2B(c, d, a, p, J2, T2z);
      double mat2 = get1B(b, p);
      double res_add = mat1 * mat2 / E_dom;
      res_E1 += res_add;
      //
      vec_d[0] += res_add;
      int Derive_phase = 1;
      double Derive_Edom = 1.0;
      for (int i = 1; i <= Derivative_Max; i++) {
        Derive_phase = Derive_phase * -1;
        Derive_Edom = Derive_Edom / E_dom;
        vec_d[i] += res_add * Derive_phase * Derive_Edom;
      }
    }
    // E2
    if ((t2z_a == t2z_p) && (j2_a == j2_p) && (l_a % 2) == (l_p % 2)) {
      double E_mid = this->hbar_omega * (e_p + e_b);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      double mat1 = get2B(c, d, p, b, J2, T2z);
      double mat2 = get1B(a, p);
      double res_add = mat1 * mat2 / E_dom;
      res_E2 += res_add;
      //
      vec_d[0] += res_add;
      int Derive_phase = 1;
      double Derive_Edom = 1.0;
      for (int i = 1; i <= Derivative_Max; i++) {
        Derive_phase = Derive_phase * -1;
        Derive_Edom = Derive_Edom / E_dom;
        vec_d[i] += res_add * Derive_phase * Derive_Edom;
      }
    }
    // E3
    if ((t2z_d == t2z_p) && (j2_d == j2_p) && (l_d % 2) == (l_p % 2)) {
      double E_mid = this->hbar_omega * (e_c + e_p);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      double mat1 = get2B(a, b, c, p, J2, T2z);
      double mat2 = get1B(d, p);
      double res_add = mat1 * mat2 / E_dom;
      res_E3 += res_add;
      //
      vec_d[0] += res_add;
      int Derive_phase = 1;
      double Derive_Edom = 1.0;
      for (int i = 1; i <= Derivative_Max; i++) {
        Derive_phase = Derive_phase * -1;
        Derive_Edom = Derive_Edom / E_dom;
        vec_d[i] += res_add * Derive_phase * Derive_Edom;
      }
    }
    // E4
    if ((t2z_c == t2z_p) && (j2_c == j2_p) && (l_c % 2) == (l_p % 2)) {
      double E_mid = this->hbar_omega * (e_p + e_d);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      double mat1 = get2B(a, b, p, d, J2, T2z);
      double mat2 = get1B(p, c);
      double res_add = mat1 * mat2 / E_dom;
      res_E4 += res_add;
      //
      vec_d[0] += res_add;
      int Derive_phase = 1;
      double Derive_Edom = 1.0;
      for (int i = 1; i <= Derivative_Max; i++) {
        Derive_phase = Derive_phase * -1;
        Derive_Edom = Derive_Edom / E_dom;
        vec_d[i] += res_add * Derive_phase * Derive_Edom;
      }
    }
  }
  res = res_E1 + res_E2 + res_E3 + res_E4;
  return res;
}

double MBPT_SEC::cal_TwF(int a, int b, int c, int d, int J2, int T2z, int Par,
                         double E_start, vector<double> &vec_d) {
  double res = 0.0;
  int e_a = SP.state_Tz[a].E;
  int l_a = SP.state_Tz[a].l;
  int j2_a = SP.state_Tz[a].j_2;
  int t2z_a = SP.state_Tz[a].tz_2;

  int e_b = SP.state_Tz[b].E;
  int l_b = SP.state_Tz[b].l;
  int j2_b = SP.state_Tz[b].j_2;
  int t2z_b = SP.state_Tz[b].tz_2;

  int e_c = SP.state_Tz[c].E;
  int l_c = SP.state_Tz[c].l;
  int j2_c = SP.state_Tz[c].j_2;
  int t2z_c = SP.state_Tz[c].tz_2;

  int e_d = SP.state_Tz[d].E;
  int l_d = SP.state_Tz[d].l;
  int j2_d = SP.state_Tz[d].j_2;
  int t2z_d = SP.state_Tz[d].tz_2;

  double res_F1 = 0.0;
  double res_F2 = 0.0;
  double res_F3 = 0.0;
  double res_F4 = 0.0;

  for (int h = 0; h <= sp_core_Tz; h++) {
    int e_h = SP.state_Tz[h].E;
    int l_h = SP.state_Tz[h].l;
    int j2_h = SP.state_Tz[h].j_2;
    int t2z_h = SP.state_Tz[h].tz_2;

    // F1
    if ((t2z_b == t2z_h) && (j2_b == j2_h) && (l_b % 2) == (l_h % 2)) {
      double E_mid = this->hbar_omega * (e_b + e_c + e_d - e_h);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      double mat1 = get2B(a, h, c, d, J2, T2z);
      double mat2 = get1B(h, b);
      double res_add = mat1 * mat2 / E_dom;
      res_F1 += res_add;
      //
      vec_d[0] += res_add;
      int Derive_phase = 1;
      double Derive_Edom = 1.0;
      for (int i = 1; i <= Derivative_Max; i++) {
        Derive_phase = Derive_phase * -1;
        Derive_Edom = Derive_Edom / E_dom;
        vec_d[i] += res_add * Derive_phase * Derive_phase;
      }
    }
    // F2
    if ((t2z_a == t2z_h) && (j2_a == j2_h) && (l_a % 2) == (l_h % 2)) {
      double E_mid = this->hbar_omega * (e_a + e_c + e_d - e_h);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      double mat1 = get2B(h, b, c, d, J2, T2z);
      double mat2 = get1B(h, a);
      double res_add = mat1 * mat2 / E_dom;
      res_F2 += res_add;
      //
      vec_d[0] += res_add;
      int Derive_phase = 1;
      double Derive_Edom = 1.0;
      for (int i = 1; i <= Derivative_Max; i++) {
        Derive_phase = Derive_phase * -1;
        Derive_Edom = Derive_Edom / E_dom;
        vec_d[i] += res_add * Derive_phase * Derive_Edom;
      }
    }
    // F3
    if ((t2z_d == t2z_h) && (j2_d == j2_h) && (l_d % 2) == (l_h % 2)) {
      double E_mid = this->hbar_omega * (e_a + e_b + e_d - e_h);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      double mat1 = get2B(c, h, a, b, J2, T2z);
      double mat2 = get1B(h, d);
      double res_add = mat1 * mat2 / E_dom;
      res_F3 += res_add;
      //
      vec_d[0] += res_add;
      int Derive_phase = 1;
      double Derive_Edom = 1.0;
      for (int i = 1; i <= Derivative_Max; i++) {
        Derive_phase = Derive_phase * -1;
        Derive_Edom = Derive_Edom / E_dom;
        vec_d[i] += res_add * Derive_phase * Derive_Edom;
      }
    }
    // F4
    if ((t2z_c == t2z_h) && (j2_c == j2_h) && (l_c % 2) == (l_h % 2)) {
      double E_mid = this->hbar_omega * (e_a + e_c + e_d - e_h);
      double E_dom = E_start - E_mid;
      if (E_dom == 0.0)
        continue;
      if (std::abs(E_dom) > this->E_cut * hbar_omega)
        continue;
      double mat1 = get2B(h, b, c, d, J2, T2z);
      double mat2 = get1B(h, a);
      double res_add = mat1 * mat2 / E_dom;
      res_F4 += res_add;
      //
      vec_d[0] += res_add;
      int Derive_phase = 1;
      double Derive_Edom = 1.0;
      for (int i = 1; i <= Derivative_Max; i++) {
        Derive_phase = Derive_phase * -1;
        Derive_Edom = Derive_Edom / E_dom;
        vec_d[i] += res_add * Derive_phase * Derive_Edom;
      }
    }
  }
  res = res_F1 + res_F2 + res_F3 + res_F4;
  return res;
}

void MBPT_SEC::print_sec_OB(std::string filename) {
  std::ofstream file;
  file.open(filename);
  LEC_CONST LECs;
  file << " ---- Parameters : ----" << endl;
  // file << "\t Lambda_Cut = " << Lambda_Cut << "  (MeV)" << endl;
  file << "\t c1 = " << LECs.LEC_C1 << "  (GeV^-1)" << endl;
  file << "\t c3 = " << LECs.LEC_C3 << "  (GeV^-1)" << endl;
  file << "\t c4 = " << LECs.LEC_C4 << "  (GeV^-1)" << endl;
  file << "\t cD = " << LECs.LEC_CD << endl;
  file << "\t cE = " << LECs.LEC_CE << endl;
  file << "\t cE = " << LECs.LEC_gA << endl;
  file << "\t Hbar-Omega = " << hbar_omega << endl;
  // file << "\t The Power of Regulator : " << reg_pow << endl;
  file << " \t +++++++++++++++++" << endl;
  file << "\t    E_cut = " << E_cut << endl;
  file << " \t +++++++++++++++++" << endl;
  file << " ----------------" << endl;

  int p_beg_Tz = sp_P_beg_Tz;
  int p_end_Tz = sp_P_end_Tz;
  int h_core_Tz = sp_core_Tz;
  int h_core = std::floor(h_core_Tz / 2.0);
  int p_beg = std::floor(p_beg_Tz / 2.0);
  int p_end = std::floor(p_end_Tz / 2.0);
  /*
    file << "  === Single-Particle Orbits === " << endl;
    file << "\t index\t n\t l\t j2\t";
    for (int i = p_end; i >= p_beg; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Model Space" << endl;
    file << "\t --- --- --- --- \t";
    for (int i = h_core; i >= 0; i--) {
      file << endl << "\t" << i + 1 << "\t";
      file << SP.state[i].n << "\t";
      file << SP.state[i].l << "\t";
      file << SP.state[i].j_2 << "\t";
    }
    file << "\t Core Space" << endl;
    file << endl;/**/

  file << "  === Single-Particle Orbits === " << endl;
  file << "\tindex\t n\t l\t 2*j\t ";
  for (int i = p_end; i >= p_beg; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[i].n << "\t";
    file << SP.state[i].l << "\t";
    file << SP.state[i].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Model Space" << endl;
  file << "\t --- --- --- --- \t";
  for (int i = h_core; i >= 0; i--) {
    file << endl << "\t" << i + 1 << "\t";
    file << SP.state[i].n << "\t";
    file << SP.state[i].l << "\t";
    file << SP.state[i].j_2 << "\t";
    // file << SP.state_Tz[i].tz_2 << "\t";
  }
  file << "\t Core Space" << endl;
  file << endl; /**/

  file << " === 5 Second Order One-Body Diagrams value === " << endl;
  file << "a\t t2z\t\t A\t\t B\t\t C\t\t D\t\t  E*2\t " << endl;
  for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
    int i = std::floor(a / 2.0) + 1;
    file << i << "\t " << SP.state_Tz[a].tz_2 << std::setw(8) << "\t "
         << std::setprecision(5) << vec_OA[a] << "\t " << std::setw(8)
         << std::setprecision(5) << vec_OB[a] << "\t " << std::setw(8)
         << std::setprecision(5) << vec_OC[a] << "\t " << std::setw(8)
         << std::setprecision(5) << vec_OD[a] << "\t " << std::setw(8)
         << std::setprecision(5) << vec_OE[a] << endl;
  }
  file << " === Normal-Ordering Single-Particle Energy === " << endl;
  file << "index\t n\t l\t 2*j\t 2*tz\t  Second_Order  " << endl;

  for (int i = p_end_Tz; i >= p_beg_Tz; i--) {
    file << std::floor(i / 2.0) + 1 << "\t";
    // file << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
    // file << std::setw(8) << std::setprecision(5) << vec_1f3_C[i] << "\t ";
    // file << std::setw(8) << std::setprecision(5) << vec_1f3_P[i] << "\t ";
    // file << std::setw(8) << std::setprecision(5) << vec_1f3_1C1P[i] <<
    // "\t\t
    // ";
    file << std::setw(8) << std::setprecision(5) << vec_1B[i] << "\t\t ";
    // file << std::setw(8) << std::setprecision(5) << vec_1f3[i] + vec_sec[i]
    //     << "\t ";
    file << endl;
  }
  file << endl;
  file.close();
} /**/

void MBPT_SEC::print_sec_OB_Luigi(std::string filename) {
  std::ofstream file;
  file.open(filename);
  LEC_CONST LECs;
  file << " ---- Parameters : ----" << endl;
  // file << "\t Lambda_Cut = " << Lambda_Cut << "  (MeV)" << endl;
  file << "\t c1 = " << LECs.LEC_C1 << "  (GeV^-1)" << endl;
  file << "\t c3 = " << LECs.LEC_C3 << "  (GeV^-1)" << endl;
  file << "\t c4 = " << LECs.LEC_C4 << "  (GeV^-1)" << endl;
  file << "\t cD = " << LECs.LEC_CD << endl;
  file << "\t cE = " << LECs.LEC_CE << endl;
  file << "\t cE = " << LECs.LEC_gA << endl;
  file << "\t Hbar-Omega = " << hbar_omega << endl;
  // file << "\t The Power of Regulator : " << reg_pow << endl;
  file << " \t +++++++++++++++++" << endl;
  file << "\t    E_cut = " << E_cut << endl;
  file << " \t +++++++++++++++++" << endl;
  file << " ----------------" << endl;

  int p_beg_Tz = sp_P_beg_Tz;
  int p_end_Tz = sp_P_end_Tz;
  int h_core_Tz = sp_core_Tz;
  int h_core = std::floor(h_core_Tz / 2.0);
  int p_beg = std::floor(p_beg_Tz / 2.0);
  int p_end = std::floor(p_end_Tz / 2.0);

  vector<State_SP> state;
  state = SP.state;
  std::stable_sort(state.begin(), state.end(), [](State_SP a, State_SP b) {
    return (100000 * a.E + 10000 * a.n - 1000 * a.l - 100 * a.j_2) <
           (100000 * b.E + 10000 * b.n - 1000 * b.l - 100 * b.j_2);
  });
  int size_sp = SP.state.size();
  int size_sp_t = state.size();
  vector<int> this2other(size_sp);
  for (int i = 0; i < size_sp; i++) {

    int n = SP.state[i].n;
    int l = SP.state[i].l;
    int j2 = SP.state[i].j_2;
    // int tz2 = SP.state_Tz[i].tz_2;
    int flag = -1;
    for (int j = 0; j < size_sp_t; j++) {
      int n_t = state[j].n;
      int l_t = state[j].l;
      int j2_t = state[j].j_2;
      //  int tz2_t = state_Tz[j].tz_2;
      if (n_t != n || l_t != l)
        continue;
      if (j2_t != j2)
        continue;
      flag = j;
      break;
    }
    if (flag < 0) {
      cout << "Wrong happened @Norm_Ord_Sec::print_2V_all" << endl;
      exit(0);
    }
    this2other[i] = flag + 1;
    // cout << "this2other[" << i << "] = " << this2other[i] << endl;
  }
  file << "n\t l\t j\t index" << endl;
  for (int i = 0; i < size_sp_t; i++) {
    int index = i + 1;
    int n = state[i].n;
    int l = state[i].l;
    int j2 = state[i].j_2;
    // int tz_2 = state_Tz[i].tz_2;
    file << n << "\t " << l << "\t " << j2 << "\t " << index << endl;
  }

  file << " === 5 Second Order One-Body Diagrams value === " << endl;
  file << "a\t t2z\t\t A\t\t B\t\t C\t\t D\t\t  E*2\t " << endl;
  for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
    int i = std::floor(a / 2.0);
    file << this2other[i] << "\t " << SP.state_Tz[a].tz_2 << "\t"
         << std::setw(8) << std::setprecision(5) << vec_OA[a] << "\t "
         << std::setw(8) << std::setprecision(5) << vec_OB[a] << "\t "
         << std::setw(8) << std::setprecision(5) << vec_OC[a] << "\t "
         << std::setw(8) << std::setprecision(5) << vec_OD[a] << "\t "
         << std::setw(8) << std::setprecision(5) << vec_OE[a] << endl;
  }
  file << " === Normal-Ordering Single-Particle Energy === " << endl;
  file << "index\t n\t l\t 2*j\t 2*tz\t  Second_Order  " << endl;

  for (int i = sp_P_beg_Tz; i <= sp_P_end_Tz; i++) {
    int a = std::floor(i / 2.0);
    file << this2other[a] << "\t";
    // file << "\t" << i + 1 << "\t";
    file << SP.state_Tz[i].n << "\t";
    file << SP.state_Tz[i].l << "\t";
    file << SP.state_Tz[i].j_2 << "\t";
    file << SP.state_Tz[i].tz_2 << "\t";
    // file << std::setw(8) << std::setprecision(5) << vec_1f3_C[i] << "\t ";
    // file << std::setw(8) << std::setprecision(5) << vec_1f3_P[i] << "\t ";
    // file << std::setw(8) << std::setprecision(5) << vec_1f3_1C1P[i] <<
    // "\t\t
    // ";
    file << std::setw(8) << std::setprecision(5) << vec_1B[i] << "\t\t ";
    // file << std::setw(8) << std::setprecision(5) << vec_1f3[i] + vec_sec[i]
    //     << "\t ";
    file << endl;
  }
  file << endl;
  file.close();
} /**/
void MBPT_SEC::print_U(string filename) {
  // string filename1 = "data/vec_U.dat";
  // string filename1="data/vec_UT.dat";
  std::ofstream out_file;
  out_file.open(filename);
  out_file << "index\t n\t l\t j2\t tz*2\t " << endl;
  vector<State_SP> state_Tz;
  state_Tz = SP.state_Tz;
  int size_sp = state_Tz.size();

  for (int i = 0; i < size_sp; i++) {
    int index = i;
    int n = state_Tz[i].n;
    int l = state_Tz[i].l;
    int j2 = state_Tz[i].j_2;
    int tz_2 = state_Tz[i].tz_2;
    out_file << index << "\t " << n << "\t " << l << "\t " << j2 << "\t "
             << tz_2 << endl;
  }
  out_file << "--++++++++++++++++++++++++++++++++++++++++++++++++++--" << endl;
  out_file << "i\t j\t val_U " << endl;
  int sp_size = SP.state_Tz.size();
  for (int i = 0; i < sp_size; i++) {
    for (int j = 0; j < sp_size; j++) {
      int flag = 0;
      double val_U = norm_ct.get1B(i, j, flag);
      out_file << i << "\t " << j << "\t " << val_U << "\t " << flag << endl;
    }
  }
  out_file.close();
}
double MBPT_SEC::cal_X0(int j2a, int j2b, int J2ab, int j2c, int j2d, int J2cd,
                        int J2ac, int J2bd) {
  double res = 0.0;
  if (J2ab != J2cd || J2ac != J2bd) {
    cout << "Wrong @ MBPT_SEC::cal_X0" << endl;
    cout << "J2ab != J2cd || J2ac != J2bd" << endl;
    exit(0);
  }
  double res_6 = this->sixJ_x.quick(j2a, j2b, J2ab, j2d, j2c, J2ac);
  int phase_t = (j2b + J2ab + j2c + J2ac) / 2;
  int phase = phase_cal(phase_t);
  res = phase * Sqrt_vec[J2ab + 1] * Sqrt_vec[J2ac + 1] * res_6;
  return res;
}
double MBPT_SEC::cal_CRC_ME_I(int a, int b, int c, int d, int J2, int T2z) {
  // cross-coupled ME <ab|V|cd>
  // 1: Jc->Ja  Jd->Jb
  double res = 0.0;
  double pre_fac = 1 / Sqrt_vec[J2 + 1];
  int j2_a = SP.state_Tz[a].j_2;
  int j2_b = SP.state_Tz[b].j_2;
  int j2_c = SP.state_Tz[c].j_2;
  int j2_d = SP.state_Tz[d].j_2;

  int t2z_a = SP.state_Tz[a].tz_2;
  int t2z_b = SP.state_Tz[b].tz_2;
  int t2z_c = SP.state_Tz[c].tz_2;
  int t2z_d = SP.state_Tz[d].tz_2;

  if (T2z != (t2z_a + t2z_b)) {
    cout << "Wrong @  MBPT_SEC::cal_CRC_ME_I A" << endl;
    exit(0);
  }

  if (T2z != (t2z_c + t2z_d)) {
    cout << "Wrong @  MBPT_SEC::cal_CRC_ME_I B" << endl;
    exit(0);
  }

  int J2min_ab = std::abs(j2_a - j2_b);
  int J2min_cd = std::abs(j2_c - j2_d);
  int J2min = std::max(J2min_ab, J2min_cd);

  int J2max_ab = std::abs(j2_a + j2_b);
  int J2max_cd = std::abs(j2_c + j2_d);
  int J2max = std::min(J2max_ab, J2max_cd);
  if (J2min > J2max) {
    // return res;
    cout << "Wrong @  MBPT_SEC::cal_CRC_ME_I" << endl;
    exit(0);
  }
  double res_t = 0.0;
  // cout << "\t --- J2:" << J2 << "  T2z:" << T2z << "\t a:" << a << "  b:" <<
  // b
  //      << "  c:" << c << "  d:" << d << endl;
  for (int J2p = J2min; J2p <= J2max; J2p = J2p + 2) {
    double C2 = cal_X0(j2_c, j2_a, J2, j2_d, j2_b, J2, J2p, J2p);
    double mat = get2B(a, b, c, d, J2p, T2z);
    double C1 = Sqrt_vec[J2p + 1];
    res_t += C1 * C2 * mat;

    // cout << "\t\t +++ J2p:" << J2p << "\t C=" << C << "\t mat=" << mat <<
    // endl;
  }
  res = pre_fac * res_t;
  return res;
}

void MBPT_SEC::print_sec_TwB_deriva() {
  // cout << " this->hbar_omega = " << this->hbar_omega << endl;

  // cout << "T2z\t J2\t a\t b\t c\t d\t\t 1st\t\t ph\t\t A\t\t D\t\t E\t\t
  // F\t\t "
  //         "B\t\t C"
  //      << endl;

  // int T2z = 0; // <pn|V|pn>
  for (int T2z = -2; T2z <= 2; T2z = T2z + 2) {
    for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
      int j2_a = SP.state_Tz[a].j_2;
      int t2z_a = SP.state_Tz[a].tz_2;
      for (int b = a; b <= sp_P_end_Tz; b++) {
        int j2_b = SP.state_Tz[b].j_2;
        int t2z_b = SP.state_Tz[b].tz_2;
        int J2_min = std::abs(j2_a - j2_b);
        int J2_max = j2_a + j2_b;
        if (T2z != t2z_a + t2z_b)
          continue;
        for (int c = sp_P_beg_Tz; c <= sp_P_end_Tz; c++) {
          int j2_c = SP.state_Tz[c].j_2;
          int t2z_c = SP.state_Tz[c].tz_2;
          for (int d = c; d <= sp_P_end_Tz; d++) {
            int j2_d = SP.state_Tz[d].j_2;
            int t2z_d = SP.state_Tz[d].tz_2;

            // if (a != 3 || b != 2 || c != 3 || d != 2)
            //   continue;
            if (T2z != t2z_c + t2z_d)
              continue;
            for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
              if (Tri(J2, j2_c, j2_d) < 0)
                continue;
              int ch = DP_Model.FindC_Tz(J2, T2z);
              int index_F = DP_Model.FindDP_Tz(ch, a, b);
              int index_I = DP_Model.FindDP_Tz(ch, c, d);
              double val_1st = get2B_normalized_Tz(a, b, c, d, J2, T2z);
              vector<vector<double>> vec_d(8);
              for (int i = 0; i < 8; i++) {
                vec_d[i].resize(Derivative_Max + 1);
                for (int j = 0; j <= Derivative_Max; j++)
                  vec_d[i][j] = 0.0;
              }
              vector<double> vec_d_2Bf1B(Derivative_Max + 1, 0.0);
              double val_2Bf1B =
                  this->cal_2Bfr1B_normalized_Tz(a, b, c, d, J2, vec_d_2Bf1B);
              //
              double val_TA = vec_TwA_d[ch][index_F][index_I][0];
              double val_TB = vec_TwB_d[ch][index_F][index_I][0];
              double val_TC = vec_TwC_d[ch][index_F][index_I][0];
              double val_T_ph = val_TB + val_TC;
              double val_TD = vec_TwD_d[ch][index_F][index_I][0];
              double val_TE = vec_TwE_d[ch][index_F][index_I][0];
              double val_TF = vec_TwF_d[ch][index_F][index_I][0];
              cout << " --- T2z=" << T2z << "\t J2=" << J2 << "\t a:" << a
                   << "\t b:" << b << "\t c:" << c << "\t d:" << d << "\t --- "
                   << endl;
              cout << "v_bare : " << val_1st << endl;
              cout << "\t Q      : ";
              for (int i = 0; i <= Derivative_Max; i++)
                cout << std::setw(8) << std::setprecision(5)
                     << Qm_vec[ch][index_F][index_I][i] << "   ";
              cout << endl;
              cout << "\t S 2Bf1B: ";
              for (int i = 0; i <= Derivative_Max; i++)
                cout << std::setw(8) << std::setprecision(5) << vec_d_2Bf1B[i]
                     << "   ";
              cout << endl;
              cout << "\t hf+U h : ";
              for (int i = 0; i <= Derivative_Max; i++)
                cout << std::setw(8) << std::setprecision(5)
                     << vec_TwF_d[ch][index_F][index_I][i] << "   ";
              cout << endl;
              cout << "\t hf+U p : ";
              for (int i = 0; i <= Derivative_Max; i++)
                cout << std::setw(8) << std::setprecision(5)
                     << vec_TwE_d[ch][index_F][index_I][i] << "   ";
              cout << endl;
              cout << "\tladder  : ";
              for (int i = 0; i <= Derivative_Max; i++)
                cout << std::setw(8) << std::setprecision(5)
                     << vec_TwA_d[ch][index_F][index_I][i] << "   ";
              cout << endl;
              cout << "\t   bk D : ";
              for (int i = 0; i <= Derivative_Max; i++)
                cout << std::setw(8) << std::setprecision(5)
                     << vec_TwD_d[ch][index_F][index_I][i] << "   ";
              cout << endl;
              cout << "\t ph B+C : ";
              for (int i = 0; i <= Derivative_Max; i++)
                cout << std::setw(8) << std::setprecision(5)
                     << vec_TwB_d[ch][index_F][index_I][i] +
                            vec_TwC_d[ch][index_F][index_I][i]
                     << "   ";
              cout << endl;
              cout << "\t  Sum  : ";
              for (int i = 0; i <= Derivative_Max; i++)
                cout << std::setw(8) << std::setprecision(5)
                     << vec_TwA_d[ch][index_F][index_I][i] +
                            vec_TwB_d[ch][index_F][index_I][i] +
                            vec_TwC_d[ch][index_F][index_I][i] +
                            vec_TwD_d[ch][index_F][index_I][i] +
                            vec_TwE_d[ch][index_F][index_I][i] +
                            vec_TwF_d[ch][index_F][index_I][i]
                     << "   ";
              cout << endl;
              // cout << "\t ph C : ";
              // for (int i = 0; i <= Derivative_Max; i++)
              //   cout << std::setw(8) << std::setprecision(5)
              //        << vec_TwC_d[ch][index_F][index_I][i] << "   ";
              // cout << endl;
            }
          }
        }
      }
    }
  } /**/
  // int a = 3;
  // int h = 0;
  // int c = 3;
  // int p = 8;
  // int J2 = 4;
  // T2z = 0;
  // double cmat = cal_CRC_ME_I(a, h, c, p, J2, T2z);
  // cout << " ++ J2 = " << J2 << "  T2z = " << T2z << "\t a:" << a << "  h:" <<
  // h
  //      << "  c:" << c << "  p:" << p << "\t cmat=" << cmat << endl;
  // J2 = 4;
  // double mat = get2B(a, h, c, p, J2, T2z);
  // cout << " ++ J2 = " << J2 << "  T2z = " << T2z << "\t a:" << a << "  h:" <<
  // h
  //      << "  c:" << c << "  p:" << p << "\t mat=" << mat << endl;
}
void MBPT_SEC::print_sec_TwB() {
  cout << " this->hbar_omega = " << this->hbar_omega << endl;

  cout << "T2z\t J2\t a\t b\t c\t d\t\t 1st\t\t 2st\t\t ph\t\t A\t\t D\t\t "
          "E\t\t F\t\t "
          "B\t\t C"
       << endl;

  // int T2z = 0; // <pn|V|pn>
  for (int T2z = -2; T2z <= 2; T2z = T2z + 2) {
    for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
      int j2_a = SP.state_Tz[a].j_2;
      int t2z_a = SP.state_Tz[a].tz_2;
      for (int b = a; b <= sp_P_end_Tz; b++) {
        int j2_b = SP.state_Tz[b].j_2;
        int t2z_b = SP.state_Tz[b].tz_2;
        int J2_min = std::abs(j2_a - j2_b);
        int J2_max = j2_a + j2_b;
        if (T2z != t2z_a + t2z_b)
          continue;
        for (int c = sp_P_beg_Tz; c <= sp_P_end_Tz; c++) {
          int j2_c = SP.state_Tz[c].j_2;
          int t2z_c = SP.state_Tz[c].tz_2;
          for (int d = c; d <= sp_P_end_Tz; d++) {
            int j2_d = SP.state_Tz[d].j_2;
            int t2z_d = SP.state_Tz[d].tz_2;

            // if (a != 3 || b != 2 || c != 3 || d != 2)
            //   continue;
            if (T2z != t2z_c + t2z_d)
              continue;
            for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
              if (Tri(J2, j2_c, j2_d) < 0)
                continue;
              int ch = DP_Model.FindC_Tz(J2, T2z);
              int index_F = DP_Model.FindDP_Tz(ch, a, b);
              int index_I = DP_Model.FindDP_Tz(ch, c, d);
              double val_1st = get2B(a, b, c, d, J2, T2z);
              vector<vector<double>> vec_d(8);
              for (int i = 0; i < 8; i++) {
                vec_d[i].resize(Derivative_Max + 1);
                for (int j = 0; j <= Derivative_Max; j++)
                  vec_d[i][j] = 0.0;
              }
              //
              double val_TA = vec_TwA_d[ch][index_F][index_I][0];
              double val_TB = vec_TwB_d[ch][index_F][index_I][0];
              double val_TC = vec_TwC_d[ch][index_F][index_I][0];
              double val_T_ph = val_TB + val_TC;
              double val_TD = vec_TwD_d[ch][index_F][index_I][0];
              double val_TE = vec_TwE_d[ch][index_F][index_I][0];
              double val_TF = vec_TwF_d[ch][index_F][index_I][0];
              double val_2st = val_TA + val_T_ph + val_TD + val_TE + val_TF;
              cout << T2z << "\t " << J2 << "\t " << a << "\t " << b << "\t "
                   << c << "\t " << d << "\t\t ";
              cout << val_1st << "\t ";
              cout << val_2st << "\t ";
              cout << val_T_ph << "\t ";
              cout << val_TA << "\t ";
              cout << val_TD << "\t ";
              cout << val_TE << "\t ";
              cout << val_TF << "\t ";
              cout << val_TB << "\t ";
              cout << val_TC << endl;
            }
          }
        }
      }
    }
  } /**/
  // int a = 3;
  // int h = 0;
  // int c = 3;
  // int p = 8;
  // int J2 = 4;
  // T2z = 0;
  // double cmat = cal_CRC_ME_I(a, h, c, p, J2, T2z);
  // cout << " ++ J2 = " << J2 << "  T2z = " << T2z << "\t a:" << a << "  h:" <<
  // h
  //      << "  c:" << c << "  p:" << p << "\t cmat=" << cmat << endl;
  // J2 = 4;
  // double mat = get2B(a, h, c, p, J2, T2z);
  // cout << " ++ J2 = " << J2 << "  T2z = " << T2z << "\t a:" << a << "  h:" <<
  // h
  //      << "  c:" << c << "  p:" << p << "\t mat=" << mat << endl;
}
// -------------------------------------------------------------------------- //
void EFF_KK::init() {
  this->Derivative_Max = MB_sec.Derivative_Max;

  this->sp_core = MB_sec.sp_core;
  this->sp_core_Tz = MB_sec.sp_core_Tz;
  this->sp_P_beg = MB_sec.sp_P_beg;
  this->sp_P_beg_Tz = MB_sec.sp_P_beg_Tz;
  this->sp_P_end = MB_sec.sp_P_end;
  this->sp_P_end_Tz = MB_sec.sp_P_end_Tz;
  this->sp_Q_beg = MB_sec.sp_Q_beg;
  this->sp_Q_beg_Tz = MB_sec.sp_Q_beg_Tz;
  this->sp_Q_end = MB_sec.sp_Q_end;
  this->sp_Q_end_Tz = MB_sec.sp_Q_end_Tz;
  Frac_vec.resize(Derivative_Max + 1);
  Frac_vec[0] = 1;
  for (int i = 1; i <= Derivative_Max; i++) {
    Frac_vec[i] = Frac_vec[i - 1] * i;
  }
  this->build_info_fold();
}

void EFF_KK::build_info_fold() {
  int size = Derivative_Max + 1;
  info_fold_vec.resize(size);

  vector<vector<vector<int>>> vec_x(size);
  for (int i = 0; i < size; i++) {
    int tot_num = i;
    vector<int> vec_t(size);
    for (int k0 = tot_num; k0 >= 0; k0--) {
      vec_t[0] = k0;
      for (int k1 = k0; k1 >= 0; k1--) {
        if (k0 + k1 > tot_num)
          continue;
        vec_t[1] = k1;
        for (int k2 = k1; k2 >= 0; k2--) {
          if (k0 + k1 + k2 > tot_num)
            continue;
          vec_t[2] = k2;
          for (int k3 = k2; k3 >= 0; k3--) {
            if (k0 + k1 + k2 + k3 > tot_num)
              continue;
            vec_t[3] = k3;
            for (int k4 = k3; k4 >= 0; k4--) {
              if (k0 + k1 + k2 + k3 + k4 > tot_num)
                continue;
              vec_t[4] = k4;
              for (int k5 = k4; k5 >= 0; k5--) {
                if (k0 + k1 + k2 + k3 + k4 + k5 > tot_num)
                  continue;
                vec_t[5] = k5;
                for (int k6 = k5; k6 >= 0; k6--) {
                  if (k0 + k1 + k2 + k3 + k4 + k5 + k6 > tot_num)
                    continue;
                  vec_t[6] = k6;
                  for (int k7 = k6; k7 >= 0; k7--) {
                    if (k0 + k1 + k2 + k3 + k4 + k5 + k6 + k7 > tot_num)
                      continue;
                    vec_t[7] = k7;
                    for (int k8 = k7; k8 >= 0; k8--) {
                      if (k0 + k1 + k2 + k3 + k4 + k5 + k6 + k7 + k8 > tot_num)
                        continue;
                      vec_t[8] = k8;
                      for (int k9 = k8; k9 >= 0; k9--) {
                        if (k0 + k1 + k2 + k3 + k4 + k5 + k6 + k7 + k8 + k9 >
                            tot_num)
                          continue;
                        vec_t[9] = k9;
                        for (int k10 = k9; k10 >= 0; k10--) {
                          if (k0 + k1 + k2 + k3 + k4 + k5 + k6 + k7 + k8 + k9 +
                                  k10 >
                              tot_num)
                            continue;
                          vec_t[10] = k10;
                          if (k0 + k1 + k2 + k3 + k4 + k5 + k6 + k7 + k8 + k9 +
                                  k10 ==
                              tot_num)
                            info_fold_vec[i].push_back(vec_t);
                        }
                      }
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  /*
  for (int i = 0; i < info_fold_vec.size(); i++) {
    cout << "-- i : " << i<<" --- "<<endl;
    for (int k = 0; k < info_fold_vec[i].size(); k++) {
      cout<<"\t "<<k;
      for (int j = 0; j < Derivative_Max; j++)
        cout << "\t " << info_fold_vec[i][k][j];
      cout << endl;
    }
  }/**/
}
void EFF_KK::build_1B() {
  int sp_size = SP.state_Tz.size();
  vec_1B_eff.resize(sp_size);
  for (int i = 0; i < sp_size; i++)
    vec_1B_eff[i] = 0.0;
  for (int i = sp_P_beg_Tz; i <= sp_P_end_Tz; i++) {
    // for (int i = sp_P_beg_Tz; i <= sp_P_beg_Tz; i++) {
    cout << "\t --- i : " << i << " --- " << endl;
    vec_1B_eff[i] = Heff_cal_1B(i);
  }
}
void EFF_KK::build_2B() {
  int ch_size = DP_Model.channel_Tz.size();
  vec_2B_eff.resize(ch_size);
  double val_t = Heff_cal_2B(0, 0, 0);

  for (int ch = 0; ch < ch_size; ch++) {
    int size = DP_Model.state_Tz[ch].size();
    vec_2B_eff[ch].resize(size);
    for (int i = 0; i < size; i++) {
      vec_2B_eff[ch][i].resize(size);
      for (int j = 0; j < size; j++)
        vec_2B_eff[ch][i][j] = 0.0;
    }
    for (int i = 0; i < size; i++) {
      for (int j = 0; j < size; j++) {

        // double val = Heff_cal_2B(ch, i, j);
        // vec_2B_eff[ch][i][j] = val;
        // vec_2B_eff[ch][j][i] = val;
      }
    }
  }
}

double EFF_KK::Heff_cal_1B(int a) {
  double res = 0.0;
  res = H_iter(MB_sec.Sm_vec[a]);
  return res;
}

double EFF_KK::Heff_cal_2B(int ch, int index_F, int index_I) {
  if (ch < 0 || index_F < 0 || index_I < 0) {
    cout << "Wrong @ EFF_KK::Heff_cal" << endl;
    exit(0);
  }
  double res = 0.0;
  vector<double> H2B(Derivative_Max);
  vector<double> H2Bf1B(Derivative_Max);
  // Build 2-Body Hamiltonion

  H2Bf1B = MB_sec.Sm2Bf1B_vec[ch][index_F][index_I];
  H2B = MB_sec.Qm_vec[ch][index_F][index_I];

  //
  /*double res_2b = H_iter(H2B);/**/
  double res_2bf1b_it = H_iter(H2Bf1B);/**/

  // double res_2b = H_Fold(H2B);
  double res_2bf1b_fd = H_Fold(H2Bf1B);

  // res = res_2b - res_2bf1b;
  return res;
}

double EFF_KK::H_Fold(vector<double> &Qm_vec) {
  double val = 0.0;
  vector<double> F_vec(Derivative_Max + 1);
  for (int i = 0; i <= Derivative_Max; i++) {
    int size_sum = info_fold_vec[i].size();
    double F_val = 0.0;
    cout << "\t --- i : " << i << " --- " << endl;
    for (int x = 0; x < size_sum; x++) {
      cout << "\t x:" << x << " : ";
      double F_val_t = Qm_vec[0];
      for (int m = 0; m < i; m++) {
        int Ord = info_fold_vec[i][x][m];
        F_val_t = F_val_t * Qm_vec[Ord]/Frac_vec[Ord];
        cout << "\t " << Ord;
      }
      cout << endl;
      F_val += F_val_t;
    }
    F_vec[i] = F_val;
  }
  // return val;
  for (int i = 0; i <= Derivative_Max; i++) {
    val += F_vec[i];
    cout << " i : " << i << "\t val : " << val << "\t F_i :" << F_vec[i]
         << endl;
  }
  return val;
}

double EFF_KK::H_iter(vector<double> &Qm_vec) {
  double val = 0.0;
  double H_last = Qm_vec[0];
  double H_next = 0.0;
  int k = 0;
  cout << " ----- ----- ------" << endl;
  while (k < 100) {
    H_next = H_next_cal(Qm_vec, H_last);
    double diff_abs = std::abs(H_next - H_last) / std::abs(H_last);
    cout << "k : " << k << "\t val = " << H_last
         << "\t diff_abs rate = " << diff_abs * 100 << " %" << endl;
    H_last = H_next;
    k++;
    if (diff_abs < 0.0001) {
      val = H_last;
      return val;
      break;
    }
  }
  cout << "Wrong @ EFF_KK::H_iter Can't get converge result" << endl;
  exit(0);
  // return val;
}
double EFF_KK::H_next_cal(vector<double> &Qm_vec, double H_last) {
  double val = 0.0;
  double H_pow = 1.0;
  for (int i = 0; i <= Derivative_Max; i++) {
    double val_add = Qm_vec[i] * H_pow / Frac_vec[i];
    H_pow = H_pow * H_last;
    if (std::abs(val_add) < EPS)
      break;
    val += val_add;
  }
  return val;
}
void EFF_KK::print() {

  cout << "\t -------- One-Body Effective Interaction -------- " << endl;
  cout << "n\t l\t j2\t tz2\t\t 1st\t\t 2rd\t\t Veff-2nd" << endl;
  for (int i = sp_P_beg_Tz; i <= sp_P_end_Tz; i++) {
    cout << SP.state_Tz[i].n << "\t ";
    cout << SP.state_Tz[i].l << "\t ";
    cout << SP.state_Tz[i].j_2 << "\t ";
    cout << SP.state_Tz[i].tz_2 << "\t\t ";
    double v_1st = MB_sec.get1B(i, i);
    cout << v_1st << "\t\t ";
    cout << MB_sec.Sm_vec[i][0] << "\t\t ";
    cout << vec_1B_eff[i] << endl;
  }
  cout << endl;
  cout << "\t -------- Two-Body Effective Interaction -------- " << endl;
  cout << "T2z\t J2\t a\t b\t c\t d\t\t 1st\t\t 2rd\t\t Veff-2nd " << endl;

  // int T2z = 0; // <pn|V|pn>
  for (int T2z = -2; T2z <= 2; T2z = T2z + 2) {
    for (int a = sp_P_beg_Tz; a <= sp_P_end_Tz; a++) {
      int j2_a = SP.state_Tz[a].j_2;
      int t2z_a = SP.state_Tz[a].tz_2;
      for (int b = a; b <= sp_P_end_Tz; b++) {
        int j2_b = SP.state_Tz[b].j_2;
        int t2z_b = SP.state_Tz[b].tz_2;
        int J2_min = std::abs(j2_a - j2_b);
        int J2_max = j2_a + j2_b;
        if (T2z != t2z_a + t2z_b)
          continue;
        for (int c = sp_P_beg_Tz; c <= sp_P_end_Tz; c++) {
          int j2_c = SP.state_Tz[c].j_2;
          int t2z_c = SP.state_Tz[c].tz_2;
          for (int d = c; d <= sp_P_end_Tz; d++) {
            int j2_d = SP.state_Tz[d].j_2;
            int t2z_d = SP.state_Tz[d].tz_2;

            // if (a != 3 || b != 2 || c != 3 || d != 2)
            //   continue;
            if (T2z != t2z_c + t2z_d)
              continue;
            for (int J2 = J2_min; J2 <= J2_max; J2 = J2 + 2) {
              if (Tri(J2, j2_c, j2_d) < 0)
                continue;
              int ch = DP_Model.FindC_Tz(J2, T2z);
              int index_F = DP_Model.FindDP_Tz(ch, a, b);
              int index_I = DP_Model.FindDP_Tz(ch, c, d);
              double val_1st = MB_sec.get2B_normalized_Tz(a, b, c, d, J2, T2z);
              double val_veff = vec_2B_eff[ch][index_F][index_I];

              cout << T2z << "\t " << J2 << "\t " << a << "\t " << b << "\t "
                   << c << "\t " << d << "\t\t ";
              cout << val_1st << "\t ";
              cout << MB_sec.Qm_vec[ch][index_F][index_I][0] << "\t\t ";
              cout << val_veff << endl;
            }
          }
        }
      }
    }
  } /**/
  // int a = 3;
  // int h = 0;
  // int c = 3;
  // int p = 8;
  // int J2 = 4;
  // T2z = 0;
  // double cmat = cal_CRC_ME_I(a, h, c, p, J2, T2z);
  // cout << " ++ J2 = " << J2 << "  T2z = " << T2z << "\t a:" << a << "  h:" <<
  // h
  //      << "  c:" << c << "  p:" << p << "\t cmat=" << cmat << endl;
  // J2 = 4;
  // double mat = get2B(a, h, c, p, J2, T2z);
  // cout << " ++ J2 = " << J2 << "  T2z = " << T2z << "\t a:" << a << "  h:" <<
  // h
  //      << "  c:" << c << "  p:" << p << "\t mat=" << mat << endl;
}
